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Abstract. We present a new algorithm for the statistical model check¬ 
ing of Markov chains with respect to unbounded temporal properties, 
including full linear temporal logic. The main idea is that we monitor 
each simulation run on the fly, in order to detect quickly if a bottom 
strongly connected component is entered with high probability, in which 
case the simulation run can be terminated early. As a result, our simu¬ 
lation runs are often much shorter than required by termination bounds 
that are computed a priori for a desired level of confidence on a large 
state space. In comparison to previous algorithms for statistical model 
checking our method is not only faster in many cases but also requires 
less information about the system, namely, only the minimum transition 
probability that occurs in the Markov chain. In addition, our method can 
be generalised to unbounded quantitative properties such as mean-payoff 
bounds. 


1 Introduction 

Traditional numerical algorithms for the verification of Markov chains may be 
computationally intense or inapplicable, when facing a large state space or lim¬ 
ited knowledge about the chain. To this end, statistical algorithms are used as 
a powerful alternative. Statistical model checking (SMC) typically refers to ap¬ 
proaches where (i) hnite paths of the Markov chain are sampled a finite number 
of times, (ii) the property of interest is verified for each sampled path (e.g. state 
r is reached), and (iii) hypothesis testing or statistical estimation is used to in¬ 
fer conclusions (e.g. state r is reached with probability at most 0.5) and give 
statistical guarantees (e.g. the conclusion is valid with 99% conhdence). SMC 
approaches differ in (a) the class of properties they can verify (e.g. bounded or 
unbounded properties), (b) the strength of statistical guarantees they provide 
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Table 1. SMC approaches to Markov chain verification, organised by (i) the class 
of verifiable properties, and (ii) by the required information about the Markov chain, 
where pmin is the minimum transition probability, IS] is the number of states, and A is 
the second largest eigenvalue of the chain. 
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(e.g. confidence bounds, only asymptotic convergence of the method towards the 
correct value, or none), and (c) the amount of information they require about 
the Markov chain (e.g. the topology of the graph). In this paper, we provide an 
algorithm for SMC of unbounded properties, with confidence bounds, in the set¬ 
ting where only the minimum transition probability of the chain is known. Such 
an algorithm is particularly desirable in scenarios when the system is not known 
(“black box”), but also when it is too large to construct or fit into memory. 

Most of the previous efforts in SMC has focused on the analysis of properties 
with bounded horizon |27l2(ll2bn2lTn4] . For bounded properties (e.g. state r is 
reached with probability at most 0.5 in the first 1000 steps) statistical guaran¬ 
tees can be obtained in a completely black-box setting, where execution runs of 
the Markov chain can be observed, but no other information about the chain is 
available. Unbounded properties (e.g. state r is reached with probability at most 
0.5 in any number of steps) are significantly more difficult, as a stopping crite¬ 
rion is needed when generating a potentially infinite execution run, and some 
information about the Markov chain is necessary for providing statistical guar¬ 
antees (for an overview, see Table [T]). On the one hand, some approaches require 
the knowledge of the full topology in order to preprocess the Markov chain. On 
the other hand, when the topology is not accessible, there are approaches where 
the correctness of the statistics relies on information ranging from the second 
eigenvalue A of the Markov chain, to knowledge of both the number [S'] of states 
and the minimum transition probability Prnin- 

Our contribution is a new SMC algorithm for full linear temporal logic 
(LTL), as well as for unbounded quantitative properties (mean payoff), which 
provides strong guarantees in the form of confidence bounds. Our algorithm 
uses less information about the Markov chain than previous algorithms that 
provide confidence bounds for unbounded properties—we need to know only 
the minimum transition probability pmin of the chain, and not the number of 
states nor the topology. Yet, experimentally, our algorithm performs in many 
cases better than these previous approaches (see Section 5). Our main idea is 
to monitor each execution run on the fly in order to build statistical hypotheses 
about the structure of the Markov chain. In particular, if from observing the 
current prefix of an execution run we can stipulate that with high probability 
a bottom strongly connected component (BSCC) of the chain has been entered, 
then we can terminate the current execution run. The information obtained from 
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execution prefixes allows us to terminate executions as soon as the property is 
decided with the required confidence, which is usually much earlier than any 
bounds that can be computed a priori. As far as we know, this is the first SMC 
algorithm that uses information obtained from execution prefixes. 

Finding prnin is a light assumption in many realistic scenarios and often does 
not depend on the size of the chain - e.g. bounds on the rates for reaction kinetics 
in chemical reaction systems are typically known, from a Prism language model 
they can be easily inferred without constructing the respective state space. 

Example 1. Consider the property of reaching state r in the Markov chain de¬ 
picted in Figure [T] While the execution runs reaching r satisfy the property and 
can be stopped without ever entering any Vi^ the finite execution paths with¬ 
out r, such as stuttutuut, are inconclusive. In other words, observing this path 
does not rule out the existence of a transition from, e.g., u to r, which, if ex¬ 
isting, would eventually be taken with probability 1. This transition could have 
arbitrarily low probability, rendering its detection arbitrarily unlikely, yet its 
presence would change the probability of satisfying the property from 0.5 to 1. 
However, knowing that if there exists such a transition leaving the set, its transi¬ 
tion probability is at least Prnin = 0.01, we can estimate the probability that the 
system is stuck in the set {t, u} of states. Indeed, if existing, the exit transition 
was missed at least four times, no matter whether it exits t or u. Consequently, 
the probability that there is no such transition and {t, u} is a BSCC is at least 
1 - h -Pmin)^- 



tjo.oi tjo.Ol 
Fig. 1. A Markov chain. 


This means that, in order to get 99% confidence that {t, u} is a BSCC, we 
only need to see both t and u around 500 time^ on a run. This is in stark 
contrast to a priori bounds that provide the same level of confidence, such as 
the (1/pmin)''®' = 100‘^(™) runs required by [3], which is infeasible for large m. 
In contrast, our method’s performance is independent of m. A 

Monitoring execution prefixes allows us to design an SMC algorithm for com¬ 
plex unbounded properties such as full LTL. More precisely, we present a new 
SMC algorithm for LTL over Markov chains, specified as follows: 

1 - (1 - = 1 - 0.99®°° Ri 0.993 
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Input: we can sample finite runs of arbitrary length from an unknown finite- 
state discrete-time Markov chain A4 starting in the initial stat^, and we are 
given a lower bound Prnin > 0 on the transition probabilities in A4, an LTL 
formula (p, a threshold probability p, an indifference region e > 0, and two 
error bounds a,/3 > oH 

Output: if P[Af \= p] >P + £, return YES with probability at least 1 — a, and 
if F[Ai \= p] < P — return NO with probability at least 1 — /3. 

In addition, we present the first SMC algorithm for computing the mean payoff 
of Markov chains whose states are labelled with rewards. 

Related work. To the best of our knowledge, we present the first SMC al¬ 
gorithm that provides confidence bounds for unbounded qualitative properties 
with access to only the minimum probability of the chain Pm\n, and the first SMC 
algorithm for quantitative properties. For completeness, we survey briefly other 
related SMC approaches. SMC of unbounded properties, usually “unbounded 
until” properties, was first considered in and the first approach was proposed 
in m, but observed incorrect in [S]. Notably, in [5S] two approaches are de¬ 
scribed. The first approach proposes to terminate sampled paths at every step 
with some probability pterm • In order to guarantee the asymptotic convergence 
of this method, the second eigenvalue A of the chain must be computed, which is 
as hard as the verification problem itself. It should be noted that their method 
provides only asymptotic guarantees as the width of the confidence interval con¬ 
verges to zero. The correctness of [15] relies on the knowledge of the second eigen¬ 
value A, too. The second approach of (25] requires the knowledge of the chain’s 
topology, which is used to transform the chain so that all potentially infinite 
paths are eliminated. In |5|, a similar transformation is performed, again requir¬ 
ing knowledge of the topology. The (pre)processing of the state space required 
by the topology-aware methods, as well as by traditional numerical methods for 
Markov chain analysis, is a major practical hurdle for large (or unknown) state 
spaces. In [3] a priori bounds for the length of execution runs are calculated from 
the minimum transition probability and the number of states. However, without 
taking execution information into account, these bounds are exponential in the 
number of states and highly impractical, as illustrated in the example above. 
Another approach, limited to ergodic Markov chains, is taken in m, based on 
coupling methods. There are also extensions of SMC to timed systems |0|. Our 
approach is also related to (TUT], where the product of a non-deterministic sys¬ 
tem and Biichi automaton is explored for accepting lassos. We are not aware 
of any method for detecting BSCCs by observing a single run, employing no 
directed search of the state space. 

Experimental evaluation. Our idea of inferring the structure of the Markov 
chain on the fly, while generating execution runs, allows for their early termina¬ 
tion. In Section 5 we will see that for many chains arising in practice, such as 

We have a black-box system in the sense of m, different from e.g. HZ] or m, where 
simulations can be run from any state. 

® Except for the transition probability bound Pmin, all inputs are standard, as used in 
literature, e.g. HZ]. 
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the concurrent probabilistic protocols from the Prism benchmark suite [2], the 
BSCCs are reached quickly and, even more importantly, can be small even for 
very large systems. Consequently, many execution runs can be stopped quickly. 
Moreover, since the number of execution runs necessary for a required precision 
and confidence is independent of the size of the state space, it needs not be very 
large even for highly confident results (a good analogy is that of the opinion 
polls: the precision and confidence of opinion polls is regulated by the sample 
size and is independent of the size of the population). It is therefore not surpris¬ 
ing that, experimentally, in most cases from the benchmark suite, our method 
outperforms previous methods (often even the numerical methods) despite re¬ 
quiring much less knowledge of the Markov chain, and despite providing strong 
guarantees in the form of confidence bounds. In Section [5J we also provide theo¬ 
retical bounds on the running time of our algorithm for classes of Markov chains 
on which it performs particularly well. 


2 Preliminaries 

Definition 1 (Markov chain). A Markov chain (MC) is a tuple At = {S, P, p), 
where 

— S is a finite set of states, 

— P : S' X S' —>■ [0,1] is the transition probability matrix, such that for every 
s € S it holds J2s'eS 

— pL is a probability distribution over S. 

We let pmin := min({P(s, s') > 0 | s,s' G S}) denote the smallest positive 
transition probability in Ad. A run of Ad is an inhnite sequence p = SqSi • • ■ of 
states, such that for all i > 0, P(si, s^+i) > 0; we let p\i] denote the state s^. A 
path in Ad is a finite prefix of a run of Ad. We denote the empty path by A and 
concatenation of paths tti and 7r2 by tti . 7r2 . Each path tt in Ad determines the 
set of runs Cone(7r) consisting of all runs that start with tt . To Ad we assign the 
probability space (Runs, A, P), where Runs is the set of all runs in Ad, T is the cr- 
algebra generated by all Cone(7r), and P is the unique probability measure such 
that P[Cone(soSi • • • Sfe)] = At(so) • ni=i where the empty product 

equals 1. The respective expected value of a random variable / : Runs —R is 

A non-empty set (7 C S' of states is strongly connected (SC) if for every 
s, s' G C there is a path from s to s'. A set of states (7 C S is a bottom strongly 
connected component (BSCC) of Ad, if it is a maximal SC, and for each s G C 
and s' G S \ C we have P(s, s') = 0. The sets of all SCs and BSCCs in Ad are 
denoted by SC and BSCC, respectively. Note that with probability 1, the set of 
states that appear infinitely many times on a run forms a BSCC. From now on, 
we use the standard notions of SC and BSCC for directed graphs as well. 
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3 Solution for reachability 


A fundamental problem in Markov chain verification is computing the probability 
that a certain set of goal states is reached. For the rest of the paper, let At = 
{S, P, /r) be a Markov chain and G C S he the set of the goal states in At. We 
let OG := {p € Runs | 3* > 0 : p[i\ G G} denote the event that “eventually a 
state in G is reached.” The event OG is measurable and its probability P[0G] 
can be computed numerically or estimated using statistical algorithms. Since no 
bound on the number of steps for reaching G is given, the major difficulty for any 
statistical approach is to decide how long each sampled path should be. We can 
stop extending the path either when we reach G, or when no more new states 
can be reached anyways. The latter happens if and only if we are in a BSCC 
and we have seen all of its states. 

In this section, we first show how to monitor each simulation run on the fly, 
in order to detect quickly if a BSCC has been entered with high probability. 
Then, we show how to use hypothesis testing in order to estimate P[0G]. 

3.1 BSCC detection 

We start with an example illustrating how to measure probability of reaching a 
BSCC from one path observation. 

Example 2. Recall Example [T] and Figure [TJ Now, consider an execution path 
stuttutu. Intuitively, does {t, m} look as a good “candidate” for being a BSCC 
of Ad? We visited both t and u three times; we have taken a transition from 
each t and u at least twice without leaving {t,u}. By the same reasoning as in 
Example [1] we could have missed some outgoing transition with probability at 
most (1 — Prnin)^- The structure of the system that can be deduced from this 
path is in Figure [5] and is correct with probability at least 1 — (1—A 

Now we formalise our intuition. Given a finite or infinite sequence p = 
sqSi ■ • ■, the support of p is the set p = {sq, si, • ■ •}■ Further, the graph of p 
is given by vertices p and edges {(si, Si+i) | i = 0,1,...}. 



Fig. 2. A graph of a path stuttutu. 


Definition 2 (Candidate). If a path tt has a suffix k such that 11, is a BSCC 
of the graph of n, we call k the candidate of tt. Moreover, for k G N, we call 
it a fc-candidate (of tt) if each s G k has at least k occurrences in k and the 
last element of k has at least k + 1 occurrences. A fc-candidate of a run p is a 
k-eandidate of some prefix of p. 
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Note that for each path there is at most one candidate. Therefore, we write K{'k) 
to denote the candidate of tt if there is one, and K{'k) = T, otherwise. Observe 
that each K{tt) _L is a SC in M.. 

Example 3. Consider a path tt = stuttutu, then K{tt) = {t, m}. Observe that {t} 
is not a candidate as it is not maximal. Further, K{7t) is a 2-candidate (and as 
such also a 1-candidate), but not a 3-candidate. Intuitively, the reason is that 
we only took a transition from u (to the candidate) twice, cf. Example [H A 

Intuitively, the higher the k the more it looks as if the /c-candidate is indeed a 
BSCC. Denoting by Candk{K) the random predicate of K being a fc-candidate 
on a run, the probability of “unluckily” detecting any specific non-BSCC set of 
states K as a, fc-candidate, can be bounded as follows. 

Lemma 1. For every K C S such that K ^ BSCC, and every s G K, k € N, 
¥[Candk{K) \ Os] < {1 - . 

Proof. Since K is not a BSCC, there is a state t G K with a transition to t' ^ K. 
The set of states K is a fc-candidate of a run, only if t is visited at least fc times 
by the path and was never followed by t' (indeed, even if t is the last state in the 
path, by definition of a fc-candidate, there are also at least fc previous occurrences 
of t in the path). Further, since the transition from t to t' has probability at least 
Prniri) the probability of not taking the transition fc times is at most (1 — Pmin)*- 

□ 

Example 4- We illustrate how candidates “evolve over time” along a run. Con¬ 
sider a run p = sqSoSiSo • • • of the Markov chain in Figure [H The empty and 
one-letter prefix do not have the candidate defined, sqSo has a candidate {sq}) 
then again iC(soSoSi) = T, and iC(soSoSiSo) = {sojSi}. One can observe that 
subsequent candidates are either disjoint or contain some of the previous candi¬ 
dates. Consequently, there are at most 2[S'! — 1 candidates on every run, which 
is in our setting an unknown bound. A 



Fig. 3. A family (for n € N) of Markov chains with large eigenvalues. 


While we have bounded the probability of detecting any specific non-BSCC 
set K as a fc-candidate, we need to bound the overall error for detecting a 
candidate that is not a BSCC. Since there can be many false candidates on a 
run before the real BSCC (e.g. Figure[3]), we need to bound the error of reporting 
any of them. 
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In the following, we first formalise the process of discovering candidates along 
the run. Second, we bound the error that any of the non-BSCC candidates 
becomes a Ic-candidate. Third, we bound the overall error of not detecting the 
real BSCC by increasing k every time a different candidate is found. 

We start with discovering the sequence of candidates on a run. For a run 
p = SqSi • ■ •, consider the sequence of random variables defined by K{sq ... Sj) 
for j > 0 , and let (iFi)i>i be the subsequence without undefined elements 
and with no repetition of consecutive elements. For example, for a run g = 
S 0 SiSiSiS 0 'SiS 2 S 2 • • •, we have Ki = {si}, K 2 = {so, si}, K 3 = {S 2 }, etc. Let Kj 
be the last element of this sequence, called the final candidate. Additionally, we 
define := Kj for all i > j. We describe the lifetime of a candidate. Given a 
non-final Ki, we write p = atfibijidiSi so that cxiUKi = 0 , fibi'yi = Ki, di ^ Ki, 
and K^aifi) 7 ^ Ki, K{ai(iibi) = Ki. Intuitively, we start exploring Ki in ff, Ki 
becomes a candidate in bi, the birthday of the ith candidate; it remains to be 
a candidate until di, the death of the ith candidate. For example, for the run 
g = S 0 S 1 S 1 S 1 S 0 S 1 S 2 S 2 ■ • ■ and i = 1, ai = sq, fi = si, bi = si, 71 = si, di = sq, 
= SlS 2 S 2 ^?[ 8 ]p[ 9 ] • • •. Note that the final candidate is almost surely a BSCC 
of Ad and would thus have 77 infinite. 

Now, we proceed to bounding errors for each candidate. Since there is an un¬ 
known number of candidates on a run, we will need a slightly stronger definition. 
First, observe that Candk{Ki) iff Ki is a fc-candidate of (dibi'yi. We say Ki is a 
strong k-candidate, written SCandk{Ki), if it is a fc-candidate of biji. Intuitively, 
it becomes a fc-candidate even not counting the discovery phase. As a result, 
even if we already assume there exists an ith candidate, its strong /c-candidacy 
gives the guarantees on being a BSCC as above in Lemma [T] 

Lemma 2. For every i, fc G N, we have 

¥[SCandkiK,) \ K^ ^ BSCC] < (1 -pr,in)'=. 


Proof. 


¥[SCandk{K,) \ K, ^ BSCC] 
1 


¥[K, ^ BSCC] 


¥[K, i BSCC] 


^ ¥[Ki = C,bi=s\- ¥[SCandk{C) \K, = C,bi = s] 


cesc\BSCC 

sec 


^Ki = C,b,=s]- F[Candk{C) \ Os] 


cesc\BSCC 

sec 


< 


(by Markov property) 

F\K- ^^BSCCI ^ = C,bi=s]-{l- pmin)'' (by Lemma[I|) 

^ ^ ‘ cesc\BSCC 

see 

= (1 - Prnin)'' (by P[K, e sc, b, e iC,] = 1 ) 


□ 






Algorithm 1 ReachedBSCC 
Input: path tt = sqSi • • • Pmin, <5 € (0,1] 

Output: Yes iff K{n) € BSCC 
C -s- _L, i ^ 0 
for j = 0 to n do 

if K{so • ■ • Sj) 7 ^ _L and K{so ■ ■ ■ Sj) ^ C then 

A'(so---Sj) 

i •<— i + 1 

i., , _ i—loKa 

^ -log(l-p™„) 

if i > 1 and SCANDfc. (A(7r), tt) then return Yes 
else return No 


Since the number of candidates can only be bounded with some knowledge 
of the state space, e.g. its size, we assume no bounds and provide a method to 
bound the error even for an unbounded number of candidates on a run. 

Lemma 3. For let Err be the set of runs such that for some t G N, 

we have SCandk^{Ki) despite Ki ^ BSCC. Then 

OO 

¥[£rr] < ^(1 • 


Proof. 


P[£’rr] = P U (^SCandk, {Ki) BSCC ! 

OO 

< ¥[SCandki H ^ BSCC] (by the union bound) 

OO 

= ^¥[SCandk,{Ki) \ Ki i BSCC] • F[Ki ef BSCC] 

OO 

<^F[SCandk,{Ki) \ Ki (f BSCC] 

(by Lemma [H]) 


2=1 




□ 

In Algorithm [1] we present a procedure for deciding whether a BSCC inferred 
from a path tt is indeed a BSCC with confidence greater than 1 — d. We use no¬ 
tation SCANDfe^(iC, tt) to denote the function deciding whether iC is a strong ki- 
candidate on tt. The overall error bound is obtained by setting ki = _ ) ■ 
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Theorem 1. For every 5>Q, Algorithm [7] is correct with error probability at 
most S. 

Proof. Since M is finite, the Algorithm [T] terminates almost surely. The proba¬ 
bility to return an incorrect result can be bounded by returning incorrect result 
for one of the non-final candidates, which by Lemma [3] is as follows: 

OO OO OO OO 

^(1 - = ^(1 - ^ S/T = S. 

i—1 i—l i—1 i—1 

□ 

We have shown how to detect a BSCC of a single path with desired confidence. 
In Algorithm [TJ we show how to use our BSCC detection method to decide 
whether a given path reaches the set G with confidence 1 — 5. The function 
NextState(7r) randomly picks a state according to p if the path is empty (tt = A); 
otherwise, if £ is the last state of tt, it randomly chooses its successor according 
to P(£, •). The algorithm returns Yes when tt reaches a state in G, and No when 
for some i, the ith candidate is a strong /c^-candidate. In the latter case, with 
probability at least 1 — 5, tt has reached a BSCC not containing G. Hence, with 
probability at most 5, the algorithm returns No for a path that could reach a 
goal. 

3.2 Hypothesis testing on a Bernoulli variable observed with 
bounded error 

In the following, we show how to estimate the probability of reaching a set of 
goal states, by combining the BSCC detection and hypothesis testing. More 
specifically, we sample many paths of a Markov chain, decide for each whether it 
reaches the goal states (Algorithmic), and then use hypothesis testing to estimate 
the event probability. The hypothesis testing is adapted to the fact that testing 
reachability on a single path may report false negatives. 

Let be a Bernoulli random variable, such that = 1 if and only if 
SiNGLEPATHREACH(G,Prnin, 5) = Yes, describing the outcome of Algorithmic 
The following theorem establishes that estimates P[0G] with a bias bounded 
by 5. 


Algorithm 2 SinglePathReach 
Input: goal states G of M, pmin,5 G (0,1] 

Output: Yes iff a run reaches G 
TT ■<— A 

repeat 

s <r- NextState(7r) 

TT t— TT . s 

if s G G then return Yes > We have provably reached G 

until REACHEDBSCC(7r,pmin, 5) 

return No > By Theorem (U V[K{'k) G BSCC] >1 — 5 
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Theorem 2. For every S > 0, we have P[0G] — 6 < < P[0G]. 

Proof. Since the event ()G is necessary for = 1, we have P[{)G | = 1] = 1. 

Therefore, P[X^ = 1] = P[0G, = 1] < P[0G], hence the upper bound. As for 

the lower bound, again P[A^ = 1] = P[0G, X^ = 1] = P[0G] — P[0G, X^ = 0] > 
P[0G] — d, where the last inequality follows by Theorem [1] and the definition of 
BSCC. □ 

In order to conclude on the value P[0G], the standard statistical model 
checking approach via hypothesis testing m decides between the hypothesis 
Ho : P[{>G] > p + e and Hi : P[0G] < p — e, where £ is a desired indifference 
region. As we do not have precise observations on each path, we reduce this 
problem to a hypothesis testing on the variable with a narrower indifference 
region: Hq : E[X^] > p + (£ — <5) and H[ : E[X|] < p — e, for some 6 < e. 

We define the reduction simply as follows. Given a statistical test T' for 
Hq,H[ we define a test T that accepts Hq if T' accepts Hq, and Hi otherwise. 
The following lemma shows that T has the same strength as T'. 

Lemma 4. Suppose the test F' decides between Hq and H[ with strength (a,/3). 
Then the test T decides between Hq with Hi with strength (a,/3). 

Proof. Consider type I error of T. Assume that Hq holds, which means P[{>G] > 
p + e. By Theorem [5] it follows that P[X^ = 1] > P[{>G] — (5 > p + (£ — (5), thus 
Hq also holds. By assumption the test T' accepts H[ with probability at most 
a, thus, by the reduction, F also accepts Hi with probability < a. The proof 
for type II error is analogous. □ 

Lemma|4]gives us the following algorithm to decide between Hq and Hi. We 
generate samples xq,xi, ■ ■ ■ ,Xn ^ X^ from SlNGLEPATHREACH(G,pmm, (5), and 
apply a statistical test to decide between Hq and H[. Finally, we accept Hq if 
Hq was accepted by the test, and Hi otherwise. In our implementation, we used 
the sequential probability ration test (SPRT) |24I23| for hypothesis testing. 

4 Extensions 

In this section, we present how the on-the-fiy BSCC detection can be used for 
verifying LTL and quantitative properties (mean payoff). 

4.1 Linear temporal logic 

We show how our method extends to properties expressible by linear temporal 
logic (LTL) |18| and, in the same manner, to all w-regular properties. Given a 
finite set Ap of atomic propositions, a labelled Markov chain (LMC) is a tuple 
M — {S,P, p, Ap, L), where {S,P,p) is a MC and L : S ^ 2^p is a labelling 
function. The formulae of LTL are given by the following syntax: 

(fi ::= a \ -<(p \ ip /\ p \ Xip \ pAJp 
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for a € Ap. The semantics is defined with respect to a word w € The 

ith letter of w is denoted by w[j], i.e. w = ?ii[0]w[l] • • • and we write Wi for the 
suffix + 1] ■ • • . We define 


w \= a a € ie[0] 

w 1= -ip ■<=^ not w \= p 
w \= ip Aif) w \= p and w \= ip 
w 1= X(^ wi\= p 

w ^ pUip 3k € N : Wk \= ip and \/ 0 < j < k : Wj \= p 


The set {w G (2'^^’)'^ | re ^ is denoted by L(y)). 

Given a labelled Markov chain Ai and an LTL formula p, we are interested 
in the measure F[M. \= p] := P[{p £ Runs | L{p) \= tp}], where L is naturally 
extended to runs by L{p)[i] = L{p[i]) for all i. 

For every LTL formula p, one can construct a deterministic Rabin automaton 
(DRA) A = (Q, 2^^", 7, (7o, Acc) that accepts all runs that satisfy p [2]. Here Q 
is a finite set of states, 7 : Q x 2^^ —>■ Q is the transition function, Qo € Q is the 
initial state, and Acc C 2*3 x 2*3 is the acceptance condition. A word w G (2"^^)“^ 
induces an infinite sequence A(iy) = spsi • • • G (5“, such that sp = Qo ^^nd 
7(si,w[j]) = Si+i for i > 0. We write Inf(w) for the set of states that occur 
infinitely often in A(w). Word w is accepted, if there exists a pair (E, F) G Acc, 
such that £’nlnf(w) = 0 and Fnlnf(w) ^ 0. The language L(A) of A is the set 
of all words accepted by A. The following is a well known result, see e.g. [21- 

Lemma 5. For every LTL formula p, a DRA A can be effectively constructed 
such that L(A) = L((/j). 

Further, the product of a MC A4 and DRA A is the Markov chain M G 



P'{{s,q),{s',q')) = 0 otherwise, and p'(s,q) = p{s) if ^{qo,L{s)) = q and 
p'{s,q) = 0 otherwise. Note that At 0 A has the same smallest transition prob¬ 
ability Pm in as A4. 

The crux of LTL probabilistic model checking relies on the fact that the 
probability of satisfying an LTL property p in a Markov chain Ai equals the 
probability of reaching an accepting BSCC in the Markov chain Ai 0 A,p. For¬ 
mally, a BSCC C of Ai 0 A^p is accepting if for some {E, F) G Acc we have 
C n (S' X A) =0 and C fl (S x F) 7^ 0. Let AccBSCC denote the union of all 
accepting BSCCs in Ai. Then we obtain the following well-known fact |2]: 


Lemma 6. For every labelled Markov chain Ai and LTL formula p, we have 
F[Ai |=p] =P[OAccBSCC]. 
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Algorithm 3 SinglePathLTL 

Input: DRA A = (Q, 2^^, 7 , ga, Acc), pmm, S € (0,1] 

Output: Yes iff the final candidate is an accepting BSCC 
q <— fjfo, ^— A 

repeat 

s NextState(7r) 

9 ^ j{<l,L{s)) 

TT ^ TT . (s,q) 

until REACHEDBSCC(7r,prnin, 5) > P[A'(7r) G BSCC] >1 — 5 

return 3{E, F) G Acc : K{tt) Cl (5 x £) = 0 A K{ti) n (S x F) / 0 


Since the input used is a Rabin automaton, the method applies to all lo- 
regular properties. Let be a Bernoulli random variable, such that = 1 
if and only if SlNGLEPATHLTL(A<^,Prnm, = Yes. Since the BSCC must be 
reached and fully explored to classify it correctly, the error of the algorithm can 
now be both-sided. 

Theorem 3. For every 5 > 0, P[At |= < P[At |= i^] -|- <5. 

Further, like in Section [321 we can reduce the hypothesis testing problem for 
FLq : P[At \= >p + e and Hi : P[At |= i^] <p — e 
for any 5 < e to the following hypothesis testing problem on the observable X^ 


Lf' :E[A^] >p+(e-5) and H[ :E[X^^] < p - (e - 6). 


4.2 Mean payoff 

We show that our method extends also to quantitative properties, such as mean 
payoff (also called long-run average reward). Let Xi = (S', P,/i) be a Markov 
chain and r : S —^ [0,1] be a reward function. Denoting by Si the random 
variable returning the i-th state on a run, the aim is to compute 



MP := lim E 


This limit exists (see, e.g. [111), and equals X^ceBSCC PIOC] • MPc,where MPc is 
the mean payoff of runs ending in C. Note that MPc can be computed from r and 
transition probabilities in C [161 • We have already shown how our method esti¬ 
mates PIOC]. Now we show how it extends to estimating transition probabilities 
in BSCCs and thus the mean payoff. 

First, we focus on a single path tt that has reached a BSCC C = K{tt) and 
show how to estimate the transition probabilities P(s, s') for each s, s' G C. Let 
Xs^s' be the random variable denoting the event that NextState(s) = s'. Xg^s' 
is a Bernoulli variable with parameter P(s, s'), so we use the obvious estimator 
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P(s,s') = #ss'(7r)/#s(7r), where #a{T^) is the number of occurrences of a in tt. 
If TT is long enough so that #s(7’') is large enough, the estimation is guaranteed 
to have desired precision ^ with desired confidence (1 — 6 s^s')- Indeed, using 
Hdffding’s inequality, we obtain 

P[P(s,s')-P(s,s')| >e] (1) 

Hence, we can extend the path tt with candidate C until it is long enough so that 
we have a. 1 — 5c confidence that all the transition probabilities in C are in the 
^-neighbourhood of our estimates, by ensuring that s'eC^ These 

estimated transition probabilities P induce a mean payoff MP^. Moreover, MPq- 
estimates the real mean payoff MP;^. Indeed, by ra, 

/ £ 

|MPc-MPc|<C:= 1 + — -1- (2) 

\ Prr\\n / 

Note that by Taylor’s expansion, for small £, we have C ~ 2|C'|£. 


Algorithm 4 SinglePathMP 

Input: reward function r, pmin, Ci ^ (0,1], 

Output: MPc such that |MPc — MPc| < C where C is 

the BSCC of the generated 

run 


TT A 

repeat 


TT <— TT . NextState(7r) 
if A(7r) A T then 


£=P™n((l+C)'/"""^"^'-l) 

> By Equation ([2]) 

,, , ln(2|K(,r)p)-ln(«/2) 

-s ^ - 

> By Equation ([T]) 

until REACHEDBSCC(7r,prtiin, 5/2) and SCANDfc(A(7r), 
return MPx(,r) computed from P and r 

tt) 


Algorithm d] extends Algorithm d] as follows. It divides the conhdence param¬ 
eters 5 into Sbscc (used as in Algorithm to detect the BSCC) and Sc (the 
total conhdence for the estimates on transition probabilities). For simplicity, we 
set Sbscc = Sc = S/2. First, we compute the bound £ required for £-precision 
(by Eq. d]). Subsequently, we compute the required strength k of the candidate 
guaranteeing dc-coiihdence on P (from Eq. [1]). The path is prolonged until the 
candidate is strong enough; in such a case M Pc is £-approximated with 1 — Sc 
conhdence. If the candidate of the path changes, all values are computed from 
scratch for the new candidate. 

Theorem 4. For every ^ > 0, the Algorithm^ terminates correctly with proba¬ 
bility at least \ — 5. 
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Proof. From Eq. [H by the union bound, we are guaranteed that the probability 
that none of the estimates Ps,s' is outside of the ^-neighbourhood doesn’t exceed 
the sum of all respective estimation errors, that is. Sc = s'ec Next, from 
Eq. [2]and from the fact that C is subject to Theorem [T] with confidence Sbscc, 

P(|MPc(r)-MPc(r)| >C) = 

=P{C G BSCC)P(|MP(r) - MP(r)| > Cl C* e BSCC)+ 

P{C i BSCC)P(|MP(r) - MP(r)| > Cl C* ^ BSCC) 

<1 • (5c + Sbscc • 1 = ^c + Sbscc < <5- 


□ 

Let random variable denote the value SlNGLEPATHMP(r, Prnin, C; ^)- The 
following theorem establishes relation between the mean-payoff MP and the ex¬ 
pected value of ^Mp. 

Theorem 5. For every (5, C>0, MP — C — <5< E[X|^p] < MP + C + <5. 

Proof. Let us write as an expression of random variables Y, W, Z 

X^P = P(1 - IT) + 

where 1) IT is a Bernoulli random variable, such that IT = 0 iff the algorithm cor¬ 
rectly detected the BSCC and estimated transition probabilities within bounds, 
2) Y is the value computed by the algorithm if IT = 0, and the real mean payoff 
MP when IT = 1, and 3) Z is any random variable with the range [0,1]. The 
interpretation is as follows: when IT = 0 we observe the result Y, which has 
bounded error C, and when IT = 1 we observe arbitrary Z. We note that Y, W, Z 
are not necessarily independent. By Theorem 0 ] E[IT] < 6 and by linearity of 
expectation: E[X^p] = E[T] — E[TfT] -|-E[1TZ]. For the upper bound, observe 
that E[T] < MP -I- C, E[TIT] is non-negative and EltTZ] < 6 . As for the lower 
bound, note that E[T] > MP — C, E[TIT] < (5 and E[1T^] is non-negative. □ 

As a consequence of Theorem[5l if we establish that with (1 — a) confidence W^p 
belongs to the interval [a, h\, then we can conclude with (1 — a) confidence that 
M P belongs to the interval [a — C — <5, + C + ^] • Standard statistical methods 
can be applied to find the confidence bound for W^’p. 

5 Experimental evaluation 

We implemented our algorithms in the probabilistic model checker Prism |13| . 
and evaluated them on the DTMC examples from the Prism benchmark suite 
HU. The benchmarks model communication and security protocols, distributed 
algorithms, and fault-tolerant systems. To demonstrate how our method per¬ 
forms depending on the topology of Markov chains, we also performed experi¬ 
ments on the generic DTMCs shown in Figured and Figure H] as well as on two 
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CTMCs from the literature that have large BSCCs: “tandem” nni and “gridworld” 

m- 


Table 2. Experimental results for unbounded reachability. Simulation parameters: a = 
P = e = 0.01, 5 — 0.001, pterm = 0.0001. TO means time-out, and MO means memory- 
out. Our approach is denoted by SimAdaptive here. Highlights show the best result the 
among topology-agnostic methods. 


Example 

BSCC 

SimAdaptive 

SimTermination|25| 

SimAnalvsis|25) 

MC 

name 

size Pmin 

no., max. size 

time 

time 

time 

analysis 

time 

bluetooth{4) 

149K 7.8 • lO”-" 

3K. 1 

2.6s 

16.4s 

83.2s 

80.4s 

78.2s 

bluetooth(7) 

569K 7.8 • 10"^ 

5.8K, 1 

3.8s 

50.2s 

284.4s 

281.1s 

261.2s 

bluetootli(lO) >569K 7.8 • 10“^ 

>5.8K. 1 

5.0s 

109.2s 

TO 


TO 

brp(500,500) 

4.5M 0.01 

1.5K, 1 

7.6s 

13.8s 

35.6s 

30.7s 

103.0s 

brp(2K,2K) 

40M 0.01 

4.5K, 1 

20.4s 

17.2s 

824.4s 

789.9s 

TO 

brp(10K,10K) 

>40M 0.01 

>4.5K, 1 

89.2s 

15.8s 

TO 

- 

TO 

crowds(6,15) 

7.3M 0.066 

>3K, 1 

3.6s 

253.2s 

2.0s 

0.7s 

19.4s 

crowds(7,20) 

17M 0.05 

>3K, 1 

4.0s 

283.8s 

2.6s 

l.ls 

347.8s 

crowds(8,20) 

68M 0.05 

>3K, 1 

5.6s 

340.0s 

4.0s 

1.9s 

TO 

eql(15,10) 

616G 0.5 

1, 1 

16.2s 

TO 

151.8s 

145.1s 

110.4s 

eql(20,15) 

1279T 0.5 

1, 1 

28.8s 

TO 

762.6s 

745.4s 

606.6s 

eql(20,20) 

1719T 0.5 

1, 1 

31.4s 

TO 

TO 

- 

TO 

herman(17) 

129M 7.6 • lO”" 

1, 34 

23.0s 

33.6s 

21.6s 

O.ls 

1.2s 

herman(19) 

1162M 1.9 • 10“® 

1, 38 

96.8s 

134.0s 

86.2s 

0.1s 

1.2s 

herman(21) 

lOG 4.7 • 10“’’ 

1, 42 

570.0s 

TO 

505.2s 

O.ls 

1.4s 

leader(6,6) 

280K 2.1 • lO"” 

1, 1 

5.0s 

5.4s 

536.6s 

530.3s 

491.4s 

leader(6,8) 

>280K 3.8 • 10“*' 

1, 1 

23.0s 

26.0s 

MO 

- 

MO 

leader(6,ll) 

>280K 5.6 • 10“’' 

1, 1 

153.0s 

174.8s 

MO 

- 

MO 

nand{50,3) 

IIM 0.02 

51, 1 

7.0s 

231.2s 

36.2s 

31.0s 

272.0s 

nand(60,4) 

29M 0.02 

61, 1 

6.0s 

275.2s 

60.2s 

56.3s 

TO 

nand{70,5) 

67M 0.02 

71, 1 

6.8s 

370.2s 

148.2s 

144.2s 

TO 

tandem(500) 

>1.7M 2.4 • 10““ 

1, >501K 

2.4s 

6.4s 

4.6s 

3.0s 

3.4s 

tandem(lK) 

1.7M 9.9 • 10"® 

1, 501K 

2.6s 

19.2s 

17.0s 

12.7s 

13.0s 

tandeni(2K) 

>1.7M 4.9 • 10"® 

1, >501K 

3.4s 

72.4s 

62.4s 

59.8s 

59.4s 

gridworld(300) 

162M 1 ■ 10"® 

598, 89K 

8.2s 

81.6s 

MO 

- 

MO 

gridworld(400) 

384M 1 ■ 10"® 

798, 160K 

8.4s 

100.6s 

MO 

- 

MO 

gridworld(500) 

750M 1 ■ 10"® 

998, 250K 

5.8s 

109.4s 

MO 


MO 

Figiaie) 

37 0.5 

1, 1 

58.6s 

TO 

23.4s 

0.4s 

2.0s 

FigHlS) 

39 0.5 

1, 1 

TO 

TO 

74.8.0s 

1.8s 

2.0s 

Figll20) 

41 0.5 

1, 1 

TO 

TO 

513.6s 

11.3s 

2.0s 

Fig|311K,5) 

4022 0.5 

2, 5 

7.8s 

218.2s 

3.2s 

0.5s 

1.2s 

FiglllK,50) 

4202 0.5 

2, 50 

12.4s 

211.8s 

3.6s 

0.7s 

1.0s 

Fig|llK,500) 

6002 0.5 

2, 500, 

431.0s 

218.6s 

3.6s 

1.0s 

1.2s 

FiglUlOK.S) 

40K 0.5 

2, 5 

52.2s 

TO 

42.2s 

25.4s 

25.6s 

FigllllOOK.b) 

400K 0.5 

2, 5 

604.2s 

5.4s 

TO 

- 

TO 


All benchmarks are parametrised by one or more values, which influence their 
size and complexity, e.g. the number of modelled components. We have made 
minor modihcations to the benchmarks that could not be handled directly by 
the SMC component of Prism, by adding self-loops to deadlock states and fixing 
one initial state instead of multiple. 

Our tool can be downloaded at [1]. Experiments were done on a Linux 64-bit 
machine running an AMD Opteron 6134 CPU with a time limit of 15 minutes 
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and a memory limit of 5GB. To increase performance of our tool, we check 
whether a candidate has been found every 1000 steps; this optimization does 
not violate correctness of our analysis. See Appendix [B] for a discussion on this 
bound. 

Reachability. The experimental results for unbounded reachability are shown in 
Tabled] The Prism benchmarks were checked against their standard properties, 
when available. We directly compare our method to another topology-agnostic 
method of [3S] (SimTermination), where at every step the sampled path is ter¬ 
minated with probability pterm ■ The approach of [3] with a priori bounds is not 
included, since it times out even on the smallest benchmarks. In addition, we 
performed experiments on two methods that are topology-aware: sampling with 
reachability analysis of [2S] (SimAnalysis) and the numerical mo del-checking al¬ 
gorithm of Prism (MC). Appendix lAl contains detailed experimental evaluation 
of these methods. 

The table shows the size of every example, its minimum probability, the 
number of BSCCs, and the size of the largest BSCC. Column “time” reports the 
total wall time for the respective algorithm, and “analysis” shows the time for 
symbolic reachability analysis in the SimAnalysis method. Highlights show the 
best result among the topology-agnostic methods. All statistical methods were 
used with the SPRT test for choosing between the hypothesis, and their results 
were averaged over five runs. 

Finding the optimal termination probability ptem for the SimTermination 
method is a non-trivial task. If the probability is too high, the method might 
never reach the target states, thus give an incorrect result, and if the value is too 
low, then it might sample unnecessarily long traces that never reach the target. 
For instance, to ensure a correct answer on the Markov chain in Figure [31 pterm 
has to decrease exponentially with the number of states. By experimenting we 
found that the probability pterm = 0.0001 is low enough to ensure correct results. 
See Appendix |A| for experiments with other values of pterm • 

On most examples, our method scales better than the SimTermination method. 
Our method performs well even on examples with large BSCCs, such as “tan¬ 
dem” and “gridworld,” due to early termination when a goal state is reached. 
For instance, on the “gridworld” example, most BSCCs do not contain a goal 
state, thus have to be fully explored, however the probability of reaching such 
BSCC is low, and as a consequence full BSCC exploration rarely occurs. The 
SimTermination method performs well when the target states are unreachable or 
can be reached by short paths. When long paths are necessary to reach the target. 



Fig. 4. A Markov chain with two transient parts consisting of N strongly connected 
singletons, leading to BSCCs with the ring topology of M states. 
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the probability that an individual path reaches the target is small, hence many 
samples are necessary to estimate the real probability with high confidence. 

Moreover, it turns out that our method compares well even with methods 
that have access to the topology of the system. In many cases, the running 
time of the numerical algorithm MC increases dramatically with the size of the 
system, while remaining almost constant in our method. The bottleneck of the 
Sim Ana lysis algorithm is the reachability analysis of states that cannot reach the 
target, which in practice can be as difficult as numerical model checking. 

LTL and mean payoff. In the second experiment, we compared our algorithm 
for checking LTL properties and estimating the mean payoff with the numerical 
methods of PRISM; the results are shown in Table[3] We compare against PRISM, 
since we are not aware of any SMC-based or topology-agnostic approach for 
mean payoff, or full LTL. For mean payoff, we computed 95%-confidence bound 
of size 0.22 with parameters 6 = 0.011, ( = 0.08, and for LTL we used the same 
parameters as for reachability. Due to space limitations, we report results only on 
some models of each type, where either method did not time out. In general our 
method scales better when BSCCs are fairly small and are discovered quickly. 


Table 3. Experimental results for LTL and mean-payoff properties. Simulation pa¬ 
rameters for LTL: a = = e = 0.01, S = 0.001, for mean-payoff we computed 95%- 

confidence bound of size 0.22 with <5 = 0.011, ^ = 0.08. 


name 

LTL 

property SimAdaptive time 

MC time 

Mean payoff 

name SimAdaptive time 

MC time 

bluetooth(lO) 

□0 

8.0s 

TO 

bluetooth(lO) 

3.0s 

TO 

brp(10K,10K) 

on 

90.0s 

TO 

brp(10K,10K) 

6.6s 

TO 

crowds(8,20) 

on 

9.0s 

TO 

crowds{8,20) 

2.0s 

TO 

eql(20,20) 

no 

7.0s 

MO 

eql(20,20) 

2.6s 

TO 

herman(21) 

no 

TO 

2.0s 

herman{21) 

MO 

3.0s 

leader(6,5) 

no 

277.0s 

117.0s 

leader(6,6) 

48.5 

576.0 

nand(70,5) 

no 

4.0s 

TO 

nand(70,5) 

2.0s 

294.0s 

tandem(2K) 

no 

TO 

221.0s 

tandem(500) 

TO 

191.0s 

gridworld(lOO) 

no -!■ on 

TO 

110.4s 

gridworld(50) 

TO 

58.1s 

Fig|l20) 

no -!■ no 

TO 

3.4 

Fig|a20) 

TO 

1.8s , 

FiglUlOOK.S) 

no 

348.0s 

TO 

Fig|l;i00K.5) 

79.6s 

TO 

FigllIlK.SOO) 

no 

827.0s 

2.0s 

FigHlK.bOO) 

TO 

2.0s , 


6 Discussion and conclusion 

As demonstrated by the experimental results, our method is fast on systems 
that are (1) shallow, and (2) with small BSCCs. In such systems, the BSCC is 
reached quickly and the candidate is built-up quickly. Further, recall that the 
BSCC is reported when a fc-candidate is found, and that k is increased with each 
candidate along the path. Hence, when there are many strongly connected sets, 
and thus many candidates, the BSCC is detected by a fc-candidate for a large k. 
However, since k grows linearly in the number of candidates, the most important 
and limiting factor is the size of BSCCs. 
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We state the dependency on the depth of the system and BSCC sizes formally. 
We pick (5 := I and let 


sim 


- log ^ log i/ 


log 


p-e+S 

p+e—S 


log 


1—p —£+(5 
1—p+£ —(5 


and 


ki = 


i — log S 


-l0g(l -Pmin) 


denote the a priori upper bound on the number of simulations necessary for 
SPRT |24I23| and the strength of candidates as in Algorithm [5J respectively. 


Theorem 6. Let R denote the expected number of steps before reaching a BSCC 
and B the maximum size of a BSCC. Further, let T := maxcgBSCC;s,s'eC ^[time to 
reach s' from s]. In particular, T G 0{B/p^.^^). Then the expected running time 
of Algorithms'^ and\^ is at most 


0{sim ■ kn+B ■ B - T). 

Proof. We show that the expected running time of each simulation is at most 
kfi+B ■ B - T. Since the expected number of states visited is bounded hy R + B, 
the expected number of candidates on a run is less than 2(i? + B) — 1. Since ki 
grows linearly in i it is sufficient to prove that the expected time to visit each 
state of a BSCC once (when starting in BSCC) is at most B ■ T. We order the 
states of a BSCC as si,..., Sb, then the time is at most X]i=i where b < B. 
This yields the result since R G 0{kB+B ■ B ■ T). 

It remains to prove that T < B/p^-^^. Let s be a state of a BSCC of size at 
most B. Then, for any state s' from the same BSCC, the shortest path from s to 
s' has length at most B and probability at least p^-,„. Consequently, if starting 
at s, we haven’t reached s' after B steps with probability at most 1 — p^\„, and 
we are instead in some state s" ^ s', from which, again, the probability to reach 
s' within B steps at least Hence, the expected time to reach s' from s is at 
most 

CX3 

i=l 

where i indicates the number of times a sequence of B steps is observed. The 
series can be summed by differentiating a geometric series. As a result, we obtain 
a bound B/p^. □ 

Systems that have large deep BSCCs require longer time to reach for the required 
level of confidence. However, such systems are often difficult to handle also for 
other methods agnostic of the topology. For instance, correctness of |1S] on the 
example in Figure|3]relies on the termination probability pterm being at most 1—A, 
which is less than 2“” here. Larger values lead to incorrect results and smaller 
values to paths of exponential length. Nevertheless, our procedure usually runs 
faster than the bound suggest; for detailed discussion see Appendix [Cl 

Conclusion. To the best of our knowledge, we propose the first statistical 
mo del-checking method that exploits the information provided by each simula¬ 
tion run on the fly, in order to detect quickly a potential BSCC, and verify LTL 
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properties with the desired confidence. This is also the first application of SMC 
to quantitative properties such as mean payoff. We note that for our method to 
work correctly, the precise value of is not necessary, but a lower bound is 
sufficient. This lower bound can come from domain knowledge, or can be inferred 
directly from description of white-box systems, such as the Prism benchmark. 

The approach we present is not meant to replace the other methods, but 
rather to be an addition to the repertoire of available approaches. Our method 
is particularly valuable for models that have small BSCCs and huge state space, 
such as many of the Prism benchmarks. 

In future work, we plan to investigate the applicability of our method to 
Markov decision processes, and to deciding language equivalence between two 
Markov chains. 
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Appendix 

A Detailed experiments 

Table U] shows detailed experimental result for unbounded reachability. Com¬ 
pared to Table [5] we included: 1) experiments for the SimTermination method 
with two other values of pterm, 2) we report the number of sampled paths as 
“samples,” and 3) we report the average length of sampled paths as “path length.” 
Topology-agnostic methods, such as SimAdaptive and SimTermination, cannot be 
compared directly with topology-aware methods, such as SimAnalysis and MC, 
however for reader’s curiosity we highlighted in the table the best results among 
all methods. 

We observed that in the “herman” example the SMC algorithms work unusu¬ 
ally slow. This problem seems to be caused by a bug in the original sampling 
engine of Prism and it appears that all SMC algorithms suffer equally from this 
problem. 
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Table 4. Detailed experimental results for unbounded reachability. Simulation parameters: a = j3 = e = 0.01, <5 = 0.001. TO means 
a timeout or memory out, and WRONG means that the reported result was incorrect. Our approach is denoted by SimAdaptive here. 
Highlights show the best result among all methods. 


Example 

name 

time 

SimAdaptive 
samples path length 

$imTermination,pterm = 10 
time samples path length 

SimTermination 
time samples 

Pterm - 10“" 

path length 

SimTermination 
time samples 

Pterm - 10“=^ 

path length 

time 

SimAnalysis 

samples path length analysis 

MC 

time 

bluetootli(4) 

2.6s 

243 

499 

185.0s 

43764 

387 

16.4s 

3389 

484 

6.4s 

463 

495 

83.2s 

219 

502 

80.4s 

78.2s 

bluetootli(7) 

3.8s 

243 

946 

697.4s 

106732 

604 

50.2s 

6480 

897 

10.2s 

792 

931 

284.4s 

219 

937 

281.1s 

261.2s 

bluetooth{10) 

5.0s 

243 

1391 

TO 

- 

- 

109.2s 

9827 

1292 

15.0s 

932 

1380 

TO 


- 

- 

TO 

brp(500,500) 

7.6s 

230 

3999 

3.2s 

258 

963 

13.8s 

258 

9758 

107.2s 

258 

104033 

35.6s 

207 

3045 

30.7s 

103.0s 

brp(2K,2K) 

20.4s 

230 

13000 

3.4s 

258 

1029 

17.2s 

258 

9127 

115.0s 

258 

98820 

824.4s 

207 

12167 

789.9s 

TO 

brp(10K,10K) 

89.2s 

230 

61999 

3.6s 

258 

960 

15.8s 

258 

10059 

109.4s 

258 

96425 

TO 


- 

- 

TO 

crowds(6,15) 

3.6s 

395 

879 

29.2s 

7947 

878 

253.2s 

7477 

8735 

TO 

- 

- 

2.0s 

400 

85 

0.7s 

19.4s 

crowds(7,20) 

4.0s 

485 

859 

32.6s 

9378 

850 

283.8s 

8993 

8464 

TO 

- 

- 

2.6s 

473 

98 

1.1s 

347.8s 

crowds (8,20) 

5.6s 

830 

824 

38.2s 

11405 

821 

340.0s 

10574 

8132 

TO 

- 

- 

4.0s 

793 

110 

1.9s 

TO 

eql(15,10) 

16.2s 

1149 

652 

303.2s 

28259 

628 

TO 

- 

- 

TO 

- 

- 

151.8s 

1100 

201 

145.1s 

110.4s 

eql(20,15) 

28.8s 

1090 

1299 

612.8s 

44048 

723 

TO 

- 

- 

TO 

- 

- 

762.6s 

999 

347 

745.4s 

606.6s 

eql(20,20) 

31.4s 

1071 

1401 

TO 

11408 

156 

TO 

- 

- 

TO 

- 

- 

TO 


- 

- 

TO 

herman(17) 

23.0s 

243 

30 

257.6s 

2101 

30 

33.6s 

381 

32 

29.0s 

279 

31 

21.6s 

219 

30 

0.1s 

1.2s 

hcrman(19) 

96.8s 

243 

40 

TO 

- 

- 

134.0s 

355 

38 

254.4s 

279 

40 

86.2s 

219 

38 

0.1s 

1.2s 

herman(21) 

570.0s 

243 

46 

MO 

- 

- 

TO 

- 

- 

MO 

- 

- 

505.2s 

219 

48 

0.1s 

1.4s 

lGader(6,6) 

5.0s 

243 

7 

7.6s 

437 

7 

5.4s 

258 

7 

5.0s 

258 

7 

536.6s 

219 

7 

530.3s 

491.4s 

leader(6,8) 

23.0s 

243 

7 

62.4s 

560 

7 

26.0s 

279 

7 

26.2s 

258 

7 

MO 


- 

- 

MO 

leader(6,ll) 

153.0s 

243 

7 

TO 

- 

- 

174.8s 

279 

7 

776.8s 

258 

7 

MO 


- 

- 

MO 

nand(50,3) 

7.0s 

899 

1627 

570.6s 

140880 

846 

231.2s 

21829 

4632 

TO 

- 

- 

36.2s 

1002 

1400 

31.0s 

272.0s 

nand(60,4) 

6.0s 

522 

2431 

TO 

- 

- 

275.2s 

25250 

4494 

TO 

- 

- 

60.2s 

458 

2160 

56.3s 

TO 

nand(70,5) 

6.8s 

391 

3343 

TO 

- 

- 

370.2s 

30522 

4643 

TO 

- 

- 

148.2s 

308 

3080 

144.2s 

TO 

tandem(500) 

2.4s 

243 

501 

59.6s 

43156 

394 

6.4s 

3318 

489 

2.0s 

412 

500 

4.6s 

219 

501 

3.0s 

3.4s 

tandem(lK) 

2.6s 

243 

1001 

328.4s 

114288 

632 

19.2s 

6932 

954 

3.4s 

858 

995 

17.0s 

219 

1001 

12.7s 

13.0s 

tandem(2K) 

3.4s 

243 

2001 

TO 

- 

- 

72.4s 

14881 

1811 

6.6s 

1093 

1985 

62.4s 

219 

2001 

59.8s 

59.4s 

gridworld(300) 

8.2s 

1187 

453 

214.4s 

46214 

349 

81.6s 

18678 

437 

77.4s 

16663 

449 

MO 


- 

- 

MO 

gridworld(400) 

8.4s 

1047 

543 

274.8s 

53152 

399 

100.6s 

18909 

531 

93.0s 

16674 

548 

MO 


- 

- 

MO 

gridworld(500) 

5.8s 

480 

637 

277.4s 

57263 

431 

109.4s 

18025 

605 

104.4s 

15684 

627 

MO 


- 

- 

MO 

FigElie) 

58.6s 

128 

140664 

TO 

- 

- 

TO 

- 

- 

TO 

- 

- 

23.4s 

115 

141167 

0.4s 

2.0s 

Figiais) 

TO 

- 

- 

2.8s 

258 

1015 

TO 

- 

- 

TO 

- 

- 

74.8s 

115 

537062 

1.8s 

2.0s 

Fig|a20) 

TO 

- 

- 

WRONG 


- 

TO 

- 

- 

TO 

- 

- 

513.6s 

119 

2195265 

11.3s 

2.0s 

FigEllK.S) 

7.8s 

1109 

2489 

TO 

- 

- 

218.2s 

23968 

5916 

TO 

- 

- 

3.2s 

896 

1027 

0.5s 

1.2s 

Fig^lK,50) 

12.4s 

1115 

4306 

TO 

- 

- 

211.8s 

23908 

5880 

TO 

- 

- 

3.6s 

881 

1037 

0.7s 

1.0s 

Fig0tlK,5OO) 

431.0s 

1002 

177777 

TO 

- 

- 

218.6s 

23951 

5903 

TO 

- 

- 

3.6s 

968 

1042 

1.0s 

1.2s 

Fig^lOK,5) 

52.2s 

1161 

20404 

2.6s 

258 

1072 

TO 

- 

- 

TO 

- 

- 

42.2s 

1057 

10100 

25.4s 

25.6s 

Fig^l00K,5) 

604.2s 

1331 

200399 

2.6s 

258 

981 

5.4s 

258 

9939 

TO 

- 

- 

TO 


- 

- 

TO 






















B Implementation details 


In our algorithms we frequently check whether the simulated path contains a 
candidate with the required strength. To reduce the time needed for this opera¬ 
tion we use two optimization: 1) we record SCs visited on the path, 2) we check 
if a candidate has been found every Cb > 1 steps. Our data structure records 
the sequence of SCs that have been encountered on the simulated path. The 
candidate of the path is then the last SC in the sequence. We also record the 
number of times each state in the candidate has been encountered. By using 
this data structure we avoid traversing the entire path every time we check if a 
strong fc-candidate has been reached. 

To further reduce the overhead, we update our data structure every Cb steps 
(in our experiments Cb = 1000). Figure [5] shows the impact of C'b on the running 
time for two Markov chains. The optimal value of Cb varies among examples, 
however experience shows that Cb « 1000 is a reasonable choice. 



Fig. 5. Total running time and time for processing candidates for a Markov chain in 
Figure [3] depending on the check bound Cb- 
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Fig. 6. Total running time and time for processing candidates for the ’eql(20,20)’ bench¬ 
mark depending on the check bound Cb- 


C Theoretical vs. empirical running time 

In this section, we compare the theoretical upper bound on running time given 
in Theorem |6] to empirical data. We omit the number of simulation runs (term 
sim in the theorem), and report only the logarithm of average simulation length. 
Figures [HE] and El present the comparison for different topologies of Markov 
chains. In Figure [7] we present the comparison for the worst-case Markov chain, 
which requires the longest paths to discover the BSCCs as a fc-candidate. This 
Markov chain is like the one in Figure [U but where the last state has a single 
outgoing transition to the initial state. Figure [8] suggests that the theoretical 
bound can be a good predictor of running time with respect to the depth of the 
system, however. Figure E] shows that it is very conservative with respect to the 
size of BSCCs. 
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Fig. 7 . Average length of simulations for a Markov chain like in Figure [3l but where 
the last state has a single outgoing transition to the initial state. 
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transient part length 


Fig. 8. Average length of simulations for the MC in Figure [H where M = 5 and N 
varies. 
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Fig. 9. Average length of simulations for the MC in Figure [H where N = 1000 and M 
varies. 
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